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Velocity distribution functions link the micro- and macro-level theories of fluid flow through 
porous media. Here we study them for the fluid absolute velocity and its longitudinal and lateral 
components relative to the macroscopic flow direction in a model of a random porous medium. We 
claim that all distributions follow the power exponential law controlled by an exponent 7 and a 
shift parameter uo and examine how these parameters depend on the porosity. We find that 7 has 
a universal value 1/2 at the percolation threshold and grows with the porosity, but never exceeds 2. 


The physics of viscous flows through porous media 
is important in such diverse areas of technology as oil 
recovery, energy storage, and tumor treatment [H-Q- 
Such flows, however, are notorious for their complex¬ 
ity stemming both from randomness of the medium 
and complicated interactions of different fluid particles. 
Macroscopic parameters characterizing fluid transport in 
porous media, like permeability (the ability of a porous 
system to transmit fluids) depend on a multitude of 
geometry-related parameters such as porosity, granule (or 
fracture) shape and size distribution, and specific surface 
area. This dependency, however, is nonuniversal and to a 
large extent known only through phenomenology or ap¬ 
proximate theories. 

The complete information about the flow of an incom¬ 
pressible fluid in a particular porous sample is contained 
in the velocity field. While this quantity can be stud¬ 
ied both experimentally and numerically [1, @ , it is 
sample-dependent. Therefore, to get a better insight into 
the connection between the macroscopic properties of the 
flow and the irregular structure of the medium, one needs 
mathematical tools that take into account randomness 
of the porous matrix and filter out irrelevant, sample- 
dependent information contained in the full velocity field. 
One such tool is the velocity distribution function (vdf) 
[i^ , which is the probability density function of the fluid 
velocity magnitude u or its longitudinal (ul) or trans¬ 
verse (mt) components. We will use a convenience nota¬ 
tion /, /l, and /t to denote the vdfs corresponding to 
u, Ml, and ut, respectively, and /{) and fff to denote /l 
restricted to K>o and K<o, respectively. 

Unlike the famous Maxwell-Boltzmann distribution for 
the ideal gas, vdfs for a fluid flow reflect the structure 
of the medium rather than the effect of the inter-particle 
collisions. Despite this difference, the vdfs are also closely 
related to important macroscopic parameters. For ex¬ 
ample, / and /l immediately imply the value of the hy¬ 
draulic tortuosity (r) 0, a quantity that measures the 
mean elongation of fluid paths in a porous medium 
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where the integrals are taken over the volume V of the 
porous sample. Similarly, for flows obeying Darcy’s law 


111, e.g. groundwater flows, the permeability (k) can be 


related to the mean fluid velocity along the macroscopic 
flow direction 


K = ipfi 
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where ip is the porosity of the medium, /i is the dynamic 
viscosity of the fluid, and VP is the pressure gradient. 
Thus, the vdfs can serve as a link between two macro¬ 
scopic parameters, k and r. 

Several reports on /, /l, and /t for various porous sys¬ 
tems at low Reynolds number (Re <?C 1) are already avail¬ 
able. Physically, for arguments much smaller than (m) 
their form is dominated by contributions from stagnant 
zones (dead-end pores and the volumes in the proximity 
of the fluid-solid boundary) Q, whereas for arguments 
> (m) their form reflects the properties of the conducting 
backbone, or the fluid paths carrying most of the fluid 
transport. For this reason the vdf is usually investigated 
in two physically distinct regimes: small (m < (m)) and 
high (m > (m)) fluid velocities. In the former case, the 
local fluid kinetic energy at percolation follows a power 
law 0, which implies a similar, power-law behavior for 
/(m), u <C (m). Far from percolation, however, the form 
of /(m) for small u depends on the porous matrix struc¬ 
ture [l2| and appears to be nonuniversal, therefore it will 
not be considered here. 


In contrast to the case of small velocities, the available 
results suggest the existence of some universality in the 
form of the vdfs for m > (m). The findings of different 
research groups, however, appear to be inconsistent with 
each other. On the one hand, several theoretical Q, ex¬ 
perimental Q and numerical [13, [H results suggest that 
/(m) can be approximated by a Gaussian with the max¬ 
imum shifted towards the mean fluid velocity. On the 
other hand, however, several teams reported nearly ex¬ 
ponential vdfs with the maximum located at 0. This in¬ 
cludes an experimental study on /, //(, and /t [1], as well 
as experiments 0,113 and numerical simulations [l3 . Hsj] 
for f^. Moreover, a qualitative transition from an ex¬ 
ponential to a Gaussian form of / was found for various 
sphere packings [l3]; however, in each case / peaked at 
M = 0. Finally, Siena et al. 0 suggested that follows 
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the stretched exponential function 
/l (^^l/(ml)) oc (ul/ exp [-/3 (ul/(ul))^] , (3) 


where /3 ,7 are model parameters. 

Although Eq. ([3]) encompasses both the exponential 
and Gaussian distributions, it is not applicable to the 
systems with the distribution maximum shifted from 0 . 
As a consequence, it predicts that the values of 7 can be 
much larger than 2 17| : a result not corroborated by any 
other research. 

To reconcile this difficulty, we conjecture that for u > 
{u) the velocity distribution functions follow the expo¬ 
nential power distribution 


f{u) = a exp 


u — Uo 

U-w 


7 -| 
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(and similar formulas for /{J", f ^, and /t), where a > 0 
is the normalizing factor, uq > 0 determines the location 
of the distribution peak, Uw > 0 denotes the scale factor 
corresponding to the distribution width, and 0 < 7 < 2 
is the shape factor. This is the simplest distribution that 
generalizes both the normal (7 = 2) and Laplace (7 = 1) 
distributions and allows for the shift of the distribution 
maximum from 0. In particular, in contrast to (jS)), the 
prefactor to the exponential function in (HJ does not de¬ 
pend on u. We also postulate that for / and /£j" there 
exists a threshold value of the porosity, ip *, such that 


Mo = 0 for <Pc < P < <P*, 
7 = 2 for p* < p < I, 


which reduces, by 1 , the number of unknown parameters 
in o for any given p. This number can be reduced to 2 
by noticing that each vdf is normalized to 1 . 

To verify Eq. ([3]), we examined numerically an effec¬ 
tively two-dimensional model of fibrous materials with 
the porous matrix built of identical, freely overlapping 
ob ject s randomly deposited on a regular lattice of size 
L [18|. We considered two obstacle shapes, disks and 
squares, both with the hydraulic diameter a = 8 lattice 
units (I.U.). The fluid was assumed to be incompressible 
and Newtonian, driven by a bulk force (gravity) small 
enough to ensure the creeping flow (Re <C 1). 

The basic numerical method used to solve the prob¬ 
lem was the Palabos (www.palabos.org) implementa¬ 
tion of the lattice Boltzmann method (LBM) with the 
Bhatnagar-Gross-Krook approximation for collisions and 
the numerical viscosity v =1/6 [lB - [2lj | . While the LBM 
is often used for solving flows in porous media, its accu¬ 
racy decreases when the channels in the porous matrix 
are too narrow. To verify whether this effect is signifi¬ 
cant in the model of overlapping objects, we also solved 
it with the finite difference (ED) method. In this case we 
used only square obstacles arranged so that the minimum 
channel width was 4 l.u. 0 . In both cases we used the 


periodic boundary conditions along the macroscopic fluid 
flow direction. As for the transverse direction, we applied 
the no-slip boundary conditions in the LBM and periodic 
ones for the ED. To minimize the finite-size effects, the 
lattice size, L = 1000 l.u. (LBM) and 2000 l.u. (ED) was 
chosen to ensure that L/a > 100 [ 2 ^ and the results were 
compared with those obtained for the system of size L/2. 
The simulation results were averaged over 20 independent 
porous samples for porosities p = 0.99, 0.95, 0.9,... down 
to the proximity of the percolation threshold Pc ~ 0.4968 
for the squares 23| and ~ 0.40 for the disks. The fluid ve¬ 
locity was measured at the underlying lattice nodes and 
binned to create histograms. 

Representative results obtained for overlapping disks 
using the LBM are shown in Fig. [TJ All velocities in 
this figure are normalized by (m), and hence effectively 
dimensionless. All data for f,f//,fT, including those 
not shown, can be fitted well to Eq. o constrained 
by Eq. dSD- As the porosity is increased from pc to¬ 
wards 1 , a semilog plot of / changes its shape from con¬ 
vex (subexponential), through linear (exponential), con¬ 
cave (superexponential), parabolic (Gaussian centered at 
0) and shifted parabolic (Gaussian shifted towards (u)). 
The form of /{) closely follows that of /, but f// vanishes 
faster than /£), especially far from pc- The exponent 7 
corresponding to /t also turns out to be dependent, 
though its value never reaches 2 (Fig. El). This extends 
the experimental findings of Ref. [ 6 | , where 7 « 1 was re¬ 
ported for /t at a fixed p that was chosen far from both 
Pc and 1. Note, however, that the current simulations 
cannot be used to reliably estimate 7 for /t in the limit 
of ^ 1, as in this limit the number of obstacles be¬ 
comes very small and hence a large value of L is required 
to avoid finite-size effects. It is thus possible that in this 
limit 7 tends to 2. The threshold value p* « 0.85 for / 
is close to its counterpart ~ 0.87 for f//; the accuracy of 
our simulations was insufficient to tell if they are actually 
different from each other. As expected, ug turns out a 
continuous function of p, growing from 0 for p < p* to 
(m) for (^ = 1 (Fig. El inset). Similar results were ob¬ 
tained for overlapping squares (data not shown). As for 
f//, which controls the probability distribution of nega¬ 
tive values of ul, we found that its tail is well described 
by exponential power distribution with 7 = 0.5 for all 
P (Fig- El). This conclusion was drawn from simulations 
for p < 0.85, as for higher porosities the decay of f// is 
so rapid, cf. Fig.lH that no reliable fitting of the data is 
possible. The value of 7 = 0.5 is also consistent with the 
experimental results obtained recently in Ref. 0 (Fig. [3]) 
and our simulation data for overlapping disks (data not 
shown). The contribution of f// to the transport is neg¬ 
ligible compared to that of f// only for relatively high 
porosities. Gonsequently, simplified theories that assume 
f// = 02M invalid close to pc- 

Thus, the following general picture emerges. At the 
percolation threshold all four velocity distribution func- 








3 


(p=0.45 


9=0.60 


9=0.75 


9=0.95 



Ul/{u> 


Ul/<U> 


Ul/<U> 


Ul/(u> 


10 


0 


10 





FIG. 1. (Color online) The probability distributions / (top row), /l (middle), and /t (bottom) for selected porosities </? = 
0.45, 0.60, 0.75, 0.95 (in columns, from left to right) obtained for overlapping disks. Solid lines show hts to Eq. (|4|l. 
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FIG. 2. (Color online) Exponent 7 for / (x), (+), and 

/t (□) as a function of porosity for overlapping disks. The 
symbol size roughly corresponds to the uncertainty of the 
results. Inset: uo/{u) as a function of porosity for / and 
ip> ip* ^ 0.85. 


FIG. 3. (Color online) Scaling of according to Eq. (|4]) 
with uo = 0 and 7 fixed at 0.5. Open symbols show the 
results for overlapping squares at porosity ip = 0.55, 0.7, 0.85; 
filled circles represent the experimental data of Ref. [^. The 
dashed line is the scaling function g{x) = exp{—^/x). 


tions, /,/{J’j/l , and /t decay in accordance with the 
power exponential distribution o with Mo = 0 and 
7 = 1/2. As the porosity is increased, 7 remains con¬ 
stant for /l , but increases for /,/]/, and /t. At some 


threshold porosity (f* the exponent 7 determined for / or 
/{/ reaches the maximum value of 2. As the porosity is 
increased above ip *, 7 stays fixed at 2 whereas parameter 
Mo becomes (/^-dependent and grows from 0 to (m) as ip ap- 
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proaches 1. Therefore, a vdf in a random porous medium 
is either subexponential (7 < 1 ), exponential (7 = 1 ), su¬ 
perexponential (2 > 7 > 1), or normal (7 = 2). If 7 = 2, 
the corresponding vdf depends on the porosity through 
the shift parameter uq. 


The significance of Eq. o is related to several factors. 
First, it describes the statistical properties of the micro¬ 
scopic velocity field for the velocities u > (u) that have a 
major contribution to the advective transport. Second, 
it is a relatively simple formula with only two unknown 
parameters at any given porosity. Third, the dependence 
of at least one of these parameters, 7 , on the porosity also 
appears to be fairly simple: in the model studied here it 
could be roughly approximated with two straight line seg¬ 
ments (Fig. [2]). Given this simplicity, Eq. (jd]) might be 
used to link porosity with vdf-dependent macrospcopic 
parameters like permeability or tortuosity. It should also 
be useful in studies on several open issues, like the mi¬ 
croscopic foundations of permeability and hydrodynamic 


disp ersion (longitudinal and transverse) of passive solutes 


uisp ' 


27| , the physical relevance of the tortuosity [28| , and 
properties of the conducting backbone ( 2 ^ . all with im¬ 
mediate practical applications. 


However, some important questions remain open. For 
example, to what extent our hypothesis remains valid 
for porous matrices with a complex, highly correlated 
structure [ 2 ^? Another problem is whether the value 
of 7 = 2 for / and /£j" actually implies the absence of 
long-range correlations in the velocity field above ip* @ 
or perhaps these correlations do remain and require that 
the power exponential distribution be supplemented with 
some less significant terms? 


In summary, we propose that a general form of the 
velocity distribution functions in disordered porous me¬ 
dia is given by a power exponential distribution with 
the shape factor 7 and location parameter uq such that 
1/2 < 7 < 2 and either uq = 0 or 7 = 2. Moreover, 7 
has a universal, porosity-independent value 1/2 both at 
the percolation threshold and for the negative part of the 
velocity component parallel to the macroscopic fluid flow 
direction. Our findings resolve several apparently con¬ 
flicting reports on the velocity distribution functions and 
open a new perspective on the formulation of a statistical 
theory of transport in porous media. 
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